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1. INTRODUCTION 

The National Aeronautics and Space Administration and the Department of Defense 
are actively involved in the development of a validated technology data base in the areas of 
control/structures inter— action, deployment dynamics and system performance for Large 
Space Structures (LSS). In the Control System Division of the System Dynamics 
Laboratory of the NASA/MSFC, a Ground Facility (GF), in which the dynamics and 
control system concepts being considered for LSS applications can be verified, has been 
designed and built under Dr. Henry Waites’ supervision [8]. The viability and versatility 
of this MSFC LSS ground test facility was recognized by the U.S. Air Force Wright 
Aeronautical Laboratory as a site for their Vibration Control of Space Structures (VCOSS) 
testing. 

One of the important aspects of the GF is to verify the analytical model for the 
control system design. The procedure is to describe the control system mathematically as 
well as possible, then to perform tests on the control system, and finally to factor those 
results into the mathematical model. 

However, development of a "correct" mathematical model of a system is still an art. 
In constructing large order structural models, various errors, such as modelling errors, 
parameter errors, improperly modeled uncertainties, and errors due to linearization of 
non-linear effect, create a great challenging task of determining "best" models for a 
dynamic system. It is recognized that it is conceivable that better performance will be 
anticipated when uncertainties are modeled through stochastic multiplicative and additive 
noise terms. Optimal control strategies generated under all possible parameter variations 
will definitely create more robust control systems, under controllability and observability 
conditions, than those generated by the usual approaches [15]. To avoid ad hoc 
assumptions regarding "a priori" statistics, Hyland [13,14,15] used the maximum entropy 
principle to determine a priori probability assignment induced from available data. A 
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main advantage of maximum entropy approach is that it sacrifices as little near nominal 
performance as possible while securing performance insensitivity over the likely range of 
modelling errors. 

The second issue addressed in this report is the reduction of the order of a higher 
order control plant. Usually, the principle is looking for a quadratically optimal but 
fixed— order compensator for a higher order plant in order to simplify implementation. 
Amongst the methods available in the literature, we studied methods developed by Hyland 

[16] and Wilson [34] in this project report. 

In this report, we first improved the computer program for the maximum entropy 
principle adopted in Hyland’s MEOP method [14] developed in the previous report. The 
new program then was tested against the testing problems ran by A. Gruzen [9]. It 
resulted very close match. Therefore, it is safe to say the program is successful. 

The second part of this report is centered at the theme of model reduction. Two 
methods were examined: Wilson’s model reduction method [34] and Hyland’s optimal 
projection (OP) method [14]. Design a computer program for Hyland’s OP method was 
attempted. Due to the difficulty encountered at the stage where a special matrix 
factorization technique is needed in order to obtain the required projection matrix, we were 
only able to have the program successively up to finding the LQG solution but not beyond. 
Apparently, a more thorough and deeper study of the OP method is needed. 

Numerical results along with computer programs which employed ORACLS are 
given in this report. 

This report is based on the final results of a special project conducted by Wan— Sik 
Choi who was a graduate student in the Mathematics Department at the University of 
Alabama. The project was supervised by Drs. Wei Shen Hsia and Stavros Belbas. 
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2. MAXIMUM ENTROPY MODELLING 
2.1. Maximum Entropy Method 


Consider a linear system: 

X = AX + BU + Wj 
Y = CX + u> 2 

where 

X e R n , U 6 R m , Y 6 R*, A 6 R nxn , B e R nxm , C 6 R &m , 
and 

SD(cj v u 2 ) = (v r v 2 ). 

We seek to determine a dynamic compensator 

Z = A Z 4- FY 


( 1 ) 


( 2 ) 


U = - KZ 

where ZeR n A € R 113 ™, F 6 R nx ^ and K e R 1113 ™ that minimizes the Quadratic Cost 
Function: 


J = 



+ U T R 2 U) dt 


( 3 ) 


where R and R are penalty matrices. The maximum entropy [26,27] (ME) design 
approach [11,12,13,14,15] is used to minimize J in the presence of parameter uncertainties. 


2 . 2 . 


Stratonovich Correction 

The stochastic integral I 4>(x(t), t) dx(t) can be defined in two ways. 
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Ito Integral : 


N— 1 


$(x(t),t)dx(t) = l.i.m. ^ -x(t.)] 


Stratonovich Integral : 

.b 


A— »0 


N— 1 


j=l 


f <5(x(t),t)dx(t) = l.i.m. J $ 

*' a a .n u 


A— *0 


3=1 


x(t ; ) + x(t. +1 ) ] 




[ x (t j+1 ) - 


where A = max(t. +1 — T). 

To find the relationship between two integrals, consider 


N— 1 

°A= l 
3=1 


$ 


x(t.) + x(t. +1 ) 








1 l ^[{(l-0)x(tp + 0x(t j+1 )},t.](x(t. +1 )-x(t j )] 2 ,O<0<i 

2 . , d x 
3=1 


It was shown by Stratonovich that with probability 1 

lim D. =1 f L* (x,t) b(x,t)dt. 
A--»0 A 2 J d x 


Therefore, 

j a $(x(t),t)dx(t) = £ $(x(t),t)dx(t) + 1 £ — [x(t)]b[x(t),t]dt (4) 

2 d x 

' V ' ' v ' V v / 

Stratonovich Ito correction 

where * denotes the integral in the sense of Ito. 

The relationship for the stochastic differential equations is as follows. 

Ito D.E.: dx = m[x t ,t]dt + r[x t> t]dy t 
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Stratonovich D.E.: dx t = m[x t ,t]dt + - r[x t ,t] 


d r[x t ,t] 


d x 

t 

m[x ,t] + I T[x ,t] £-L- 
1 2 1 d x t 


dt + T[x t ,t]dy t 


dt 4- T[x t ,t] dy t 


correction 


Above result was shown in [30] by using (4) and also proved in [35], 


2.3. Stochastic Modelling of Errors 

In most instances, the errors are made in the modelling process and some 

parameters may vary. Therefore, the actual system would be represented by 

p 

A actual = A+ I “i« Ai 
i=l 


where 


B actual and C 


i: set of uncorrelated uncertainties 

a(t): zero— mean, unit intensity multiplicative 

white noise 

A.; Parameter error distribution matrices 
1 

actual ta ^ e a s ' m ^ ar f° rm - 


( 5 ) 


Substituting (5) into X(t) = AX(t) yields 

p 

X(t) = (A + \ a.(t)A.) X(t) ; O.D.E. 
1=1 


P 

dx t = (A dt + ^ da. t A.)X t ; Ito S.D.E 
i=l 

p 

= AX. dt + V da. A. X 
i L it i t 

i=l 


( 6 ) 
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By comparing (6) with I t D.E. and Stratonovich D.E. we obtain 

o 

p T p 

A + I Y A 2 dt + y da.. A. X t : Stratonovich D.E. 

0 L i L it i t 

i=l J i=l 

P 

■ 1 V 1 2 

Stratonovich correction for X(t) = Ax(t) is - ^ A. 

2 i=l 

B and C take similar form. 

S B 

2.4. Necessary Conditions for Optimality [10] 

Necessary conditions take the form of two Riccati equations and two Lyapunov 
equations, all coupled by the stochastic parameters. 
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V 2s = v 2 + I c i CQ + Q) C T 


i=l 


p g = b'J’p + l B^ (P + P) A. 


i=l 


Q, = Q C I + l A i (Q + Q) C I 


i=l 


a q.=\-v£ c . 


A = A - B Rl'p 

ps s s zs s 


The compensator matrices are, 

A = \ - V#. - *K*. + D B 2 s P 
F = 

K = R^P S 
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2.5. Algorithm 

Compute F ,F • generate a stabilizing gain matrix (F) for 

initializing the solution of Riccati eq. 

Solve for • Solve Riccati eqs without having parameter 

LQG, P, Q uncertainties — uncoupled eqs. 

Begin Interations 
with LQG P, Q 

Solves P — Riccati 
no P 

converges || P,| — | P. || < c ? where | • | is a Euclidean Norm. 
? 

Solves Q— Riccati 
no Q 

converges || Qj| — I Qj < e q ? 

? 


Solves P— Lyapunov 


no P 

converges || PJ - | Pj < e- 
? 


Solves Q— Lyapunov • No need to iterate Q— Lyapunov because 

parameter doesn’t include Q 

no P,Q 

converge || P. | + | Q. | - { I P._j I } 1 < e ? 

? 

Form A c , F, k • Compensator matrices 
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2.6. Solution of Riccati equation and Lyapunov equation 

As we have seen in the necessary condition of model reductions and Maximum 

Entropy Method, the necessary conditions consist of Lyapunov equations or coupled 

Riccati and Lyapunov equations. 

Therefore solution of Riccati and Lyapunov is required for the design of control 
system. A lot of algorithm [8,18,24,28,31,32] were proposed in the past. 

In this section, algorithms which empoloyed for this special project are briefly 
discussed. 

Kleinman [19] proposed an algorithm which is based on the method of successive 
substitution to solve the algebraic Riccati equation. 

Consider the linear time-invariant system. 


X(t) = AX(t) + BU(t) X(0) = X Q 
where [A,B] is completely controllable. 

The cost to be minimized is 

00 

J(X 0 ; U(-)) = J 0 [X'(t) C' CX(t) + U'(t) R U(t)j dt 
where R is positive definite and [A,C] is completely observable. Necessary conditions for 
optimality are 

U*(X(t)) = -R -1 B / K X(t) 

and 0 = KA + A'K + C'C-KBR _1 B'K 

where K is positive definite and 

J(X ;U*(-)) = min J(X ; U(-)) = X' K X. 

U(-) 

Thus for arbitrary feedback law (X(t)) , 

J(X 0 ;U L (.)) = X ' 0 V L X . 

e (A- B L)'t ( C , C + L'RL).e (A-Bl)l dt 




eigenvalues with negative real parts. 

0 = (A - BL)'V l + VJA - BL) + C'C + L'RL. 


Kleinman’s Theorem . 

Let V k , k = 0,1, - ■ ■, be the (unique) positive definite solution of the linear algebraic 
equation 

0 = A i V k + V k A k + C 'C + L i ;RL k 
where, recursively, 

L k = R -1 B' V k _j, k = 1,2,- ■ - 
A k = A - BL k 

and where L Q is chosen such that A Q = A — BL Q has eigenvalues with negative real parts. 
Then 

1) K< V k+1 < V k < ••• , k = 0,1, • • • 

2.) lim V = K 

k-tao 

Note . In this project, stabilizing matrix L Q is computed by CSTAB in ORACLS and 
Riccati equation is solved by RICNWT in ORACLS [1]. 

An algorithm for the solution of the matrix equation AX + XB = C was proposed 
by Bartels and Stewart [6]. Above equation has a unique solution if and only if 

+ A ® f 0 (i = 1,2, - ■ ■ , m; j = 1,2,- ■ • ,n) where A^ and A® are eigenvalues of A and B 
respectively [2]. The method of solution is based on the reduction of A and B to the real 
schur form, i.e., block lower (upper) triangular form. 
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Let 


AX + XB = C 

and U, V be the orthogonal matrix. 
Then 


B' = V X BV 


B = VB'V J 


B — » upper Hessenberg form — ♦ upper real Schur form; B' 
(Heusehalder’s method) (QR algorithm) 


( 7 ) 


( 8 ) 


A' = U T A U=>A = UA' U T 


A' (lower real Schur form) is obtained by reducing the 
transparse of A to upper real Schur form and transposing back. 


C'=U 1 CV=t C = UC'V 1 
Substituting (8), (9), (10) into (7) yields 

U A' U T X + X V B' V T = U C'V T 


A' U T X + U T XVB'V T = C' V T 


A' U T X V + U T XVB'=C' 
A' X' + X' B' = C' 


c: 


ii 


pi 


( 9 ) 


( 10 ) 


A ' 
A 11 

0 


x 1 
X 11 

■ • X' 1 

iq 


x 1 • ■ 

X 11 

■ x' 
iq 



. . • b : 1 

iq 

A' A' 

z 1 zz 



• 

• 

+ 

D 



B 22 

. . . B' 

. 22q 

■ • 

• 


• 

• 





0 

■ 

A ' A ' ■ 

• • A' 


x\ • 

• ■ x' 


x' • • 

• x' 



B' 

Pi P2 

pp- 


pi 

pq J 


pi 

pq 


qq 


c: 


iq 


C' 

pq- 


k-1 t - 1 

A 'kj Xjl-l X' i B V k = l,2,.,p,k = l,2,..,q (11) 
j=l i=l 
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Equation (11) can be solved successively for X'^ . Let the right side of (11) be D. 

Since the block matrices A£ k and B' u are of order at most two, we are again required to 
solve the matrix equation of the form (7). 


Writing (11) in matrix form gives 


a ll 

a ; 2 ' 


11 

X / 
12 

+ 

x / 

X 11 

x f 
X 12 


»n 

b 12 


d n 

d i 2 ’ 

/21 

a 22. 


x / 

21 

x 1 
22 


x / 
21 

x / 

22 _ 


b 2i 

CM 
^ CM 

X 


d 21 

, 

CM 

CM 

x) 


Right side 
of (11) 


a ii + b n 


a 


21 

*12 


a i2 


b 12 


b 21 


a 22 + a U ° 


a ll + b 22 
a 21 


b 21 

a i2 

a 22 + b 2: 



1 

X 

t— ‘ 'v 
1 


d n 


x * 
x 21 

= 

d 21 


x y 

12 


d 12 

. 

x / 

L 22- 


- d 22- 


( 12 ) 


X/ .is obtained from (12). Then the solution of (7) is given by X = U X' V . 
k i 


Note . In this project, Lyapunov equation is solved by BARSTW in ORACLS [1]. 


2.7. Numerical Example for Maximum Entropy Method 

The following system posed by Doyle [9] was solved by Gruzen [10]. In this project 
some problem is solved for comparison of numerical results. 


X, 1 


■ i r 


' x .' 


' 0 


' 1 ‘ 

1 

— 



+ 


U + 


• V 


0 1. 


- X 2 ■ 


1 + nb 


1 
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Y = [1 0] 



+ V 



' 1 1 
1 1 


R. 


= 1 


V 


1 


= M 


' 1 1 
1 1 


= 1 


0, fi : parameters related with the gain margin 
Parameter uncertainty distribution matrices: 


A l = 


0 0 
0 0 


B. 


0 

fi} 


Cj = [0, 0] 


Note : 0 = fi = 60 , 0.93 < 1 + Ab < 1.01 
0 < (3 < 0.2 , size 0.05 is used. 
Necessary conditions for this example are 


0 = PA + A T P - PB T r'T B T P + R 

s s s 2s s 1 

0 = A Q + QA* - QC* vj, C Q + V, + (B, P,) Q (B, P„) T 


0 = PA. + A:F + P T r; 1 

Qs Qs s 2s 


P 

s 


0 = Ap Q + QAp + Q s V^j q'J' 

s s 


A s = A, B s = B, C s = C, R 2i = R 2 + Bf (P + P) Bj, 

V 2s = V 2’ P 5 = B I P ' Q, = ■ A Q. “ A .-Q. K <V 


where 
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A = A -B 

PS S 8 



P . 
8 


The compensator matrices are, 


A = A - Q V7 1 C - B R7 1 P 

c s 2s s s 2s 

K = K p „ 



Table 1. Numerical Results *Column II: Results for this project 



1 1 


Column I is a numerical result obtained by A. Gruzen. 
Column II is a numerical result obtained by this project. 




















18 


2.8. Discussions on ME method 

As shown in table I, matrix K decreases as 0 (disturbance) in matrix increases. 
This is because K = R“* P s , R 2s = R 2 + B^ (P + P)B 1 and similarly for matrix F but F 
increases as j3 increases. 

When 0=0 (LQG case), the two results (A.G. & N.R.) are exactly same. But for 
0t 0 best results obtained for 0 = .15. Differences in numerical results between A.G. & 
N.R. are possibly occurred from the value of Ab. (In this project Ab = 0 is used, but A. 
Gruzen doesn’t show the value of Ab which he was used). 

As a whole, the results are pretty close each other. Therefore, this indirectly 
verifies that "ME FORTRAN" provides correct answers. And it supports the fact that 
ORACLS is a good design package for designing controllers. 

3. MODEL REDUCTION: WILSON’S METHOD [34] 

3.1. Problem Statement 

Given an nth — order system 

X = AX + BU 
Y= HX, 

find an rth - order reduced system 

X = A X + B U 

r r r r 

Y = H X . 

r r r 

The input vector U(t) will be taken as a white noise, i.e., 

E(t)] = 0 

E[U(t)U T (s)] = N<5(t— s). 


(13) 

(14) 

(15) 

(16) 
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The cost function to be minimized is 

J = lim E[e T (t) Qe(t)] (17) 

t — * OD 

where e is the reduction error, e = y — y r and 

Q is positive definite. Without loss of generality assume Q is m x m idenity matrix. 


Note , where A, B, H are n x n, nxp, m x n matrices, 

A , B , H are r x r, r x p, m x r matrices, 
r r r 

x, y are n x 1, m x 1 vectors, 

x , y are r x 1, m x 1 vectors, 

r r 

U is p x 1 vector. 


3.2. Necessary conditions for optimum 


A = 0.A0, 

r 1 2 

B = 0, B 

r 1 

H r =H0 2 

where ©j = — P 2 2 ^12 an< ^ ®2 — ^12 ^"22 ' 

©l 0 2 = I r 

FR + RF T + S = 0 
F T P + PF + M = 0 


3.3. 


Derivation of Necessary Conditions 
Equation (13) - (16) may be written as 


Z = FZ + GU 
where Z = I" i 



(18) 

(19) 

( 20 ) 

( 21 ) 

( 22 ) 

(23) 


( 24 ) 
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From (17) 

J = lim E[e T Q e] 

t — * 0D 

T 

= lim E[e e] since we assumed Q = I 

t — * 00 

= lim E[(Y - Y ) T (Y - Y)] 

t — ► 0D 

= lim E[(HX - HX r ) T (HX - H r X r )] 

t — 1 0D 

Now, 

(HX - H f X r ) T (HX - H r X r ) 

= X T H T HX - X T H T H X - X T H T HX + X T H T H X 

rr rr rrrr 


rprp rprp rp rp r T' r P 

= X T H T HX-X T H T HX-X 1 H 1 HX +X i H 1 HX 

rr rr rrrr 



M 


= Z T M Z. 

Thus, 

J = lim E[Z T M Z] 

t — > 0D 
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= trace (RM) 


(25) 


where R = lim E[Z(t) Z T (t)] 

t— ♦ QD 

Let, r(t) = E[Z(t) Z T (t)]. 

Then, r(t) = E[Z(t) Z T (t) + Z(t) Z T (t)] 

= E[Z(t) Z T (t)] + E[Z(t) Z T (t)]. 

Since Z T = Z t F t + U t G a , 

r(t) = E[(FZ + GU)Z T ] + E[Z(Z T F T + U T G T )] 

= FE[ZZ T ] + GE[UZ T ] 4- E[ZZ T ]F T + E[ZU T ]G T 

= F r(t) + r(t) F T + GE[UZ T ] + E[ZU T ]G T . (26) 

But, 

Z(t) = $(t,t o ) Z(t o ) + f t 4(t,A) G(A) U(A) dA 

0 

where $(t,t) is the state transition matrix. 

Thus, 


E[UZ t ] = E[U(t) Z T (t o )] $ T (t,t o ) + f t E[U(t) U t (A)] G t $ T (t,A) dA 

o 



0 
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= J t $(t,A) G(A) Ntf(A-t) dA (28) 

o 

Substituting (27) and (28) into (26) yields 

t t 

r(t) = Fr(t) + r(t) F T + J t GM (t - A) G T $ T (t,A)dA + J t $(t,A)G(A) N<5 (A - t) G T dA 

0 0 

= Fr(t) + r(t) F T + I G N G T * T (t,t) + I $(t,t) G N G T 

2 2 

= Fr(t) + r(t) F T + G N G T . 

Since R = lim r(t), FR + RF T + G N G T = 0. 

t — > 0D 

BNB'J' ‘ 

B NB T 

r r J 

(29) 

To minimize (25) subject to (29) form the 
Lagrangian 

L = tr[ARM] + (FR + RF T + S)P. 

l_L = o => AM + F T P + PF = 0 
d R 

Let A = 1. Then 

F T P + PF + M = 0 (30) 

By comparing (30) with (29) we may write 
J = trace (PS) 


Let S = G N G 1 = 


Then, 


BNB 


B NB 

r 


FR + RF 1 + S = 0 


( 31 ) 
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Let the symmetric matrices P and R be partitioned as 

p_ P 11 P 12 R _ R 11 R 12 

K 2 P 2 J ’ R 22- ' 

Differentiating J with respect to any parameter (3 , 


Ll = 2 trf — RP 


j + tr 


Lip] + trf — R 


To find A , obtain derivative of J with respect to a using (32). Then 

r r 


i = 2 trf_ RP] + trf _ 
i I d a da, 

r L r J r 


LI Pi + trfLM R 


= 2 tr 


0 0 


d A RP 

r 


where R= R n R i2 and P 

- R 1 2 R 22 - 


P P 
11 12 

P T P 
12 22 


= 2 tr 3 A d A 

1 R 1 2 1 R 22 P T 

1-1-5 a 12 . a J 12 

r o a 

r 


1 (R 12 P 11 + R 2 /[ 2 ) 


d A / r t p p p 
' n T2' 12 + 22 22 


_2tr 7) ^ R 12 P 12 + R 22 P 22 

3 a 

L r 



24 


0 => R 12 P 12 + R 22 P 22 ° 


(33) 


P 12 R 12 + P 22 R 22 ® 


P -1 P T R + R =0 
22 12 12 ^ 22 


From (29) 


(34) 


r a o ' 


o 

> 



R 11 R 12 
R 12 R 22 


+ 


R 


12 


11 

T R 
12 22 


r a t 0 


0 A 


BNB T BNB^ 

B NB T B NB T 

r r r 


= 0 


AR 
A R 


l l 
T 

r 12 


AR 

j 

F 

r fv 22 J 


12 
A R n 


+ 


R n A 

r t a t 
R 12 A 


R 12 A T 

R 22 A T 


bnb t bnb'J' 

B NB T B NB T 

r r r 


= 0 


AR 12 + R 12 A T + BNB T = ° 
A R 22 + R 22 A ! + B , NB T = 0 


(35) 


B «t B t = - p^p> 

Thus (35) becomes 

AR 12 + R 12 A T - BNBTP 12 P ^ = 0 

A r R 22 + R 22 A ! + P ; 2 P T 2 BNBTP 12 P 22 = 0 


(36) 

(37) 


Now, T n (36) + (37) gives 


P 22 P 12 AR 12 + A r R 22 + ^ P 22 P 1 2 R 12 + R 22 ^ A r 0 


-1«T 


0, by (34) 


A r — P 22 P 12 A R 12 R 22 


same as sq. (18) 
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To find B , LA = 0. 
r d b 


d J = 2 tr 

[LZ hpI 

+ tr 

'Ll P 

+ tr 

[r* m 1 

d b 

d b 


d b 


d b 

r 

r 


r 


r 


tr 

\ p d SI 

= tr 

p — 

rBNB T BNB T 1 

r 



d b 

r - 


d b 

! r 

B NB T B NB T 

L r r r J 



= tr 


■ p n 

P 12 


'0 

BN ' 

= tr 

• p i 2 BN 

P u BN + 2P 12 B , N ' 



ip t 
L 12 

P 22- 


BN 

2B N 

r J 

J 

- P 22 BN 

P ! 2 BN + 2P 22 B , N - 


= tr (P 12 BN + P^BN + 2P 22 B r N) = tr (P^BN + P^BN + 2P 22 B r N) 
-* 2P 22 BN = -2P^BN 

B r = - P 22 P T2 B = ©i B f— same as ( 19 ) 

To find H , 12 = 0. 
r d h 


'd 

M r] 

= tr 

' d 

h t h 

H T H 1 

r 

R 

d 

h 

r J 


d h 

r 

- h t h 

L r 

h t h 

r r J 



= tr 



- H 
2H 


11 
T 

12 


R 

R 


12 

22 -* 


-tr (- HR -T 2 - HR 12 + 2H r R 22^ 


= tr (-2HR 12 + 2HR 22 ) 


=> h = hr 12 r-; = H0 2 - 



same as (20) 
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From (2.21), R 12 P 12 F"22 P 22 


pip p T > 

12 12 — 22 22 ' 


Now, 


®1 ®2 P 22 P 12 R 12 R 22 


“ P 22 ^ P 22 R 22^ R 22 


= I same as (21) 


When the conditions on A f , B r and H r {(18), (19), (20)} substituted into eqns. (29) and 
(30), a set of nonlinear equations in the unknown matrices ©j and © 2 is obtained. Namely, 

R 22 0 I AT0 1 + 0 1 A0 2 R 22 + H0 2 N0 I HT = ° 


P 22 0 1 A0 2 + 0 I A X P 22 + 0 I H H0 2=° 


An explicit solution for ©j and © 2 is not apparently possible. © L and © 2 are nonunique, in 
the sense that the output of the reduced model is invariant under any nonsingular 
transformation T. 

An algorithm to solve this optimum reduced order model problem was presented by 
Mishra and Wilson [22], 
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3.4. Algorithm [22] 

Step 1: Choose the matrices Q and N 

Step 2: Choose a value for the parameter A satisfying 0 < A < 1. Normally, 

without prior knowledge choose A = 1. 

Step 3: Make initial guesses for the matrices and B r , such that the pair (A f ,B r ) 
defines a completely controllable, strictly stable system. 

Step 4: Solve the matrix equation FT + RF T + S = 0 

Step 5: Compute the matrix 0 2 = R 12 R 22 

Step 6: Set H r = H© 2 

Step 7: Solve the matrix equation F T P + PF + M = 0 

Step 8: Compute the matrix © 1 = — P^P 12 

Step 9: Set B f = ©jB 

Step 10: If B f computed in Step 9 is not the same as B f used in Step 4, then go to 
Step 4 using the B r from Step 9. Otherwise, the B f computed in Step 9 and the 
H computed in Step 6 are taken to be the optimum for the present A f matrix. 
Step 9 and the H f computed in Step 6 are taken to be the optimum for the 
present A^ matrix. 

Step 11: Compute the error function J using the present A r matrix and the 
optimum B f and defined in Step 10. 

Step 12: Designate the present A f matrix as A° ld and the present value of the error 
function as Jq. 

Step 13: Compute a new A f . 

A" ew = A Q x AQ 2 + (1 - A)A° ld 

where ©^ and © 2 were used to compute the optimum B^ and H r for A° ld . 

Step 14: If (A“ ew , B r ) is strictly stable controllable, then go to Step 15. Otherwise, 
reduce A and go to Step 13. 
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Step 15: For A" ew and the optimum B f for A° ld , use Steps 4 to 10 until the 
optimum and H r are obtained for A" ew . 

Step 16: Compute J using A" ew , B f and H r defined in Step 10. Designate the value 
of J as Jj. 

Step 17: Test 

(a) If J x < J Q : Go to Step 12 

(b) if J > J Q : Decrease A and go to Step 13 

(c) If J = J : If 0 0=1 step. The triple (A new , B , H ) used to compute 

XU 11 r r r r 

are the optimal reduced model. Otherwise decrease A and go to Step 13. 

3.5. Derivatives of Cost Function. 


J = tr(RM) (25) 

FR 4- RF T + S = 0 (29) 

J = tr(PS) (30) 

F T P + PF + M = 0 (31) 


9 J = tr 

9 r m 

+ tr 

R (9 M 

, where @ is any parameter 

d / 3 

[d P J 


d P\ 



tr 

9 R (F T P + PF) 

+ tr 

[r 5 m 1 

since M = — (F T P + PF) from (31) 


d 0 


d P\ 



2 tr 

9 r pf 

+ tr 

[r* m 1 


[« p J 


9 P\ 


( 38 ) 
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Differentiating (29) with respect to 0, 


12 R + F 12 + 12 F t + R - I— + — = 0 (39) 

50 5050 50 50 


Postmultiply (39) by P and taking the trace 


9 F RP + F - — P + - R F T P + R - F P + 12 P = 0 
50 50 50 50 50 


tr 


12 RP 
[d 0 J 

+ tr 

p 5 Rp 

1 d 0 J 

+ tr 

'5 RpTp 

[9 P J 

+ tr 

[r d fT p] 

0 0 J 

+ tr 

'5 S p 
9 0 . 



tr 


5 R 


PF 


[5 0 


tr 

5 Rrp 


[9 0 J 


2 tr 

5 R pp 

= 2 tr 

[LErpI 

+ tr 

[5 S pi 


[9 0 \ 


[» 0 J 


[0 0 J 


Substituting (40) into (38), 


11= 2 tr 

12 RP 

+ tr 

'5 S p 

+ tr 

r! 

M] 

5 0 

1* 0 \ 


9 0 . 


5 

/?J 


= 0 


( 40 ) 


1 


4. MODEL REDUCTION: HYLAND’S METHOD [16]. 
4.1. Problem Statement 
Given the system 

X = AX + BU 
Y= CX 


find a reduced — order model 
X = A X + B U 

r r r r 

Y = C X 

r r r 


(41) 

(42) 

(43) 

( 44 ) 
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which minimizes the model — reduction criterion 


J(A,B,,C) = lira E[(Y - y/r(Y - Y)]. 


(45) 


t * GO 


The input U(t) is taken to be white noise with positive - definite intensity V. 


Note . A, B, C: 


A , B , C : 

r r r 


R, V: 


x, u, y, x f , y r : 


Pi 2 )- 


nxn, nxm, hn matrices 
nxn nxm, £ x n matrices 

r r, r r 

i x t, m x m p.d. matrices 
n, m, £, n r , t dimensional vectors 
rank of matrix Z 


Assumption: A, A^ stable. 


4.2. Necessary Conditions for Optimum 


A = rAG T 

r 

(46) 

B = TB 

r 

(47) 

c = cg t 

r 

(48) 

d Q) = d?) = p( QP) = n 

(49) 

0 = AQ + QA T + BVB T - 7i BVB T 7^ 

(50) 

0 = A T P + PA + C T RC - 7^C T RC7 x 

(51) 

where G = Q-'Q^ , T = - P^, 


7 =G T r , 7 ± = I n -7- 


TG t = I 



r 


4.3. Derivation of Necessary Conditions 
Introducing the augmented system 

X = A X + BU, 

Y= C X 
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where 



J(A B C ) = lim E[(Y - Y) T R( Y - Y )] 

t— > OD 

= tr QR where R = C T RC and Q = lim E[X(t)X T (t)]. (52) 

t — * CD 

As shown in Wilson’s Method (25) ~ (29) Q is given by the unique solution of 

0 = AQ + QA T + V (53) 

where V = BVB T 

To minimize (52) subject to (53), form the 
Lagrangian L(A r , B r , C r , Q) = tr[AQ R + (A Q + Q A T + V)P] 
where A > 0 and P € ® (n + V x (n + - 

Expanding L(A r , B f , C , Q) gives 

L = tr \\(Q 1 C T RC - Q 12 C^RC - Q^C^RC^ + Q^RCy 

+ AQjP, + AQ l2 P^ + aq; 2 p 12 + aq 2 p 2 

+ Ql ATp ! + Q 12 A T P T 2 + + VT P 2 


+ BVB T P 1 + BYB^P^ + B r VB T P 12 + B^^Pj. 
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And, 


^1 ^12 

, P = 

rp p i 
1 12 

, r = c t rc = 

[C t RC - C t RC 1 

r 

■^12 ^2 - 


p T p 
L 12 2 J 


- c t rc c t rc 

l r r r J 


V = BVB T = 


BVB 


B VB 

u r 


BVB 


B VB 

r r- 


Now, 


d Q 



d L 

d L ' 

d L = 

d Qj 

5 Ql2 

d Q 

d L 

a l 


- d 

5 q 2 . 


AC T RC +A t P 1 +P l A -AC T RC r -l- A t P 12 + P^ 


- AC T RC + A T P'[_ + p'La AC t RC +A t P + P A 

r rl2 12 rr r 2 2r 



'C t RC 

- C T RC 1 

r 


■a t Pi 

_ a t P 1 
A *12 


Pi A 

p i 2 a ; 

A 



+ 



+ 




i 

k 

- H 
J* 

O 

c t rc 

r r J 


ax. 

A T p 
r 2 J 


1 

hd 
^ H 
> 

V,- 



= AR + 



0 




= AR + A T P + PA 
Thus, A T P + PA + AR = 0. 

Without loss of generality, take A = 1. 
Then A T P + PA + R = 0 

LI i=o, 

d A 

r 


14 = 2 P » Q 1 * + 2P 2« 2 

o A 

r 

Thus, P} 2 Q 12 + P 2 Q 2 = 0 =» Q^ 2 P 12 + Q 2 P 



d B 

r 


?-L = pT 0 BV + pT 0 BV + 2P B V 
d B 12 

r 


Thus, 2[P'[ 2 B + P 2 B r ]V = 0 


d C 

r 

LL=- rcq 12 - rcq 12 + 2R C[ Q 2 

0 L> 


Thus 


2R[C r Q 2 -CQ, 2 ] = 0 
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Define, 

G = and r = - P- 1 P^ . 

Then, 

pg t = - P; lp T 2 Q I2 Qf • 

But from (55), P^Q^ = - p 2 Q 2 = - P 2 Q 2 ■ 
Thus, 

pg t = -p 2 i (-p 2 qX T = >„ r 


From (56), = — P^P^B = 

From (57), C f = CQ^Q" 1 = C(Q” T Q^ 2 ) T , Q 2 is P.d. 

= C(Q-QT/ = cg t 


Expanding (53) and (54) yields 


0 = AQ X + QjA T + BVB T 

(58) 

0 = AQ 12 + QX + BVB I 

(59) 

0 = A Q, + QjA? + B VB* 

(60) 

0 = A T P, + PjA + C T RC 

(61) 

0 = A Tp i 2 + p 12 A -C T RC 

(62) 

0 = A?P 2 + P 2 A t + C>C 

(63) 


Since A r , B f and C r are independent of Qj and (58) and (61) can be ignored. 
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Define Q = = Q 12 G 

( 64 ) 

P = P l2 P X 2 = - P 12 r 

Now (64) -r T yields 

( 65 ) 

gr T = q 12 g r T = Q 12 (rc T ) T = q 12 . 

( 66 ) 

Similarily, from (65) 


p, 2 =-pg t . 

( 67 ) 

rQP T = -P; l P^Q 12 Q; I QT 2 (-P 12 Pf) 

= Q 2 


Thus, Q 2 = TQr T 

( 68 ) 

Similarily, P 2 = GPG T 

( 69 ) 

Substitute (47), (48), (66), - (69) into (59), (60), (62), (63) 


^ rn A rp rp rp rp 

0 = AQr T + Qr aJ + BVbT 1 

( 70 ) 

o = A r rQr T + rQr T A^ + rBVB T r T 

(71) 

0 = A t PG t + PG 1 A f + C i RCG 1 

( 72 ) 

0 = A T GPG T + GPG T A + gc t rcg t . 

r r 

( 73 ) 


(71) — r-(70), 


a rQr T = rAQr T 


Q, 


X 
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Thus, A = TAQ^Q ' 1 = rAG T 

TriA f r\ 1 \ / r> — 


7Q — g tq — ( q 12 Q 2 )(p 2 p 12 )(Q 12 Q 2 Q , 2 ) 


- P 2«2 




^12*^2 *^12 ^ 


Similarily, P7 = P 
Finally, G T *(70) T yields 

G T rQ T A T + G T A r rQ T + G T rBV T B T = 0 


TAG' 


7 QA T + G T rAG T rQ + 7BV B A = 0 , Q and V symmetric. 



7 [AQ + QA T + BVB T ] = 0 


Similarily, (72) -T yields 


[A t P + PA + C A RC] 7 = 0 


(74) 

(75) 


(76) 

(77) 


(76) + (76) T + (76 ) -7 

= 7AQ + 7QA T + 7BVB T + QA T 7 T + AQ7 T + BVB T 7 T + 7AQ7 + 7 QA T 7 + 
7BVB T 7 

= QA T + AQ + 7BVB T + BVB T 7 T + 7AQ7 T + 7QA T 7 T + 7AQ7 + 7 QA T 7 + 
7BVB T 7 

= AQ + QA T + 7BVB T + BVB T 7 T + 7(AQ + QA T )7 T + 7 (AQ + QA T )7 + 7BVB T 7 
= AQ + QA T + 7BVB T + BVB T 7 T - 7BVB T 7 T 
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^ a rp rp rp rp rp rp rp rp rp 

= AQ + QA t + BVB t - BVB t + 7 BVB T r + ^BVBS - 7BVB7 

a a rp rp rp rp rp rp rp rp rp rp 

= AQ + QA T + BVB T - (I n BVB T r - 7 BVB r r - IhBVBS 1 + 7 BVB A 7 ) 

= AQ + QA + BVB T - 7 x BVB T 7 ± 

which is the same as (50) 

Similarily, 

(77) + (77) T + 7^(77) = A T P + PA + C T RC - 7^C T RC7 x 
which is the same as (51). 

A computer program has been designed (appendix 3) for this algorithm. Due to the 
difficulty of finding the projection matrix r through a matrix factorization process, the 
program only run successively up to obtaining an LQG solution. Apparently, more words 
and researchs need to be done in that area. 


4.4. Algorithm ([17,7]) 

Step 1: Initialize 7^ = I . 

Step 2: Solve for Q (K) , P (K) from 

0 = (A - 7 (K) A 7 [ K) ) <^ (K) + <^ (K) (A - 7 (K) A7[ K) ) T + BVB T 

0 = (A - 7 x (K) A7 (K) ) T £ (K) + £ {K) (A - 7 | K) A7 (K) ) + C T RC 
Step 3: Balance 

*( K >& K > ( *(K)) T = ($W)- T pw^r 1 = s (K) , 

S< K > = diag (<r[ K \ ... , ^), „{*> > a™ > • • ■ > a™ > 0 
Step 4: If K > 1 check for convergence 


e 


k 


f tr(C T RCW c ) - t r ( C T RC7 (K) Q (K) (7 (K)T ) 


l 1 / 


2 


tr(C T RCW c ) 
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If |e k — e k _J < tolerance then go to step 8), else continue. 
Step 5: Select eigenprojections 

n. [4 <K )#< K >], ..., n. n 

1 m 

n.[4' K >p< K >] i *( K )E.(i'< K) r 1 . 

N 

m 

Step 6: Update 7 (K+1) = l n. [Q (K) P (K) ] 

r=l r 

Step 7: Increment K and return to Step 2. 

Step 8: Set Q = 7^Q(7^) T > P = ( 7 ^) T P 7^ 


4.5. Relationship between two methods 


Wilson’s Method Hyland’s Method 


X = AX + BU 

X = AX + BU 

Y= HX 

Y= CX 

X = A X + B U 

r r r r 

X = A X + B U 

r r r r 

Y = H X 

r r r 

Y = C X 

r r r 

J = limE[(Y-Y [ ) T (Y-Y)) 

t — * OD 

J = lim E[(Y - Y r ) T R(Y - 

t — ► OD 

A = 0 A 0„ 

r 1 2 

A = rAG T 

r 

B = 0.B 

r 1 

B = TB 

r 

H = H0 2 

C r = cg t 

r = - p i lp T 2 

®2 = R 12 R 22 

g t = q 12 q-’ 

0 l 0 2 = 

rG T = i 

n 

r 


7 C T I 

FR + RF T + S = 0 

AQ + QA T + V = 0 



Wilson’s Method 



39 


Hyland’s Method 

A T P + PA + R = 0 
bvb t bvb t 

V = r 

B VB T B VB T 

L r r r 

c t rc - c t rc 

R = 

- c'J’rc c^rc 

AQ + QA T + BVB T - 7 i BVB T 7 '[ = 0 
A T P + PA + C T RC - 7^C T RC7 i = 0 
where 7=1—7 

'j. n 

i) necessary to solve n x n Lyapunov 
equation [50, 51] which is independent 
of A , B , and C 

r r r 

ii) Need eigenprojections to form 

N 

m 

7= l n.[QP 

i =1 

iii) Need (G, M, T) - factorization 
of Q P to determine G and T. 
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APPENDIX 1 


Numerical Results of M.E. Method 



FORTVS ME ( GS OPT ( 2 ) 

• VS FORTRAN COMPILER ENTERED. 22:45:04 

**MAIN** END OF COMPILATION 1 ****** 

**SUB1** END OF COMPILATION 2 ****** 

**SUB5** END OF COMPILATION 3 ****** 

**SUB8** END OF COMPILATION 4 ****** 

**SUB9** END OF COMPILATION 5 ****** 

**SUB12** END OF COMPILATION 6 ****** 

**SUB13** END OF COMPILATION 7 ****** 
VS FORTRAN COMPILER EXITED. 22:45:07 


GLOBAL TXTLIB VFORTLIB CMSLIB FORTUTIL 

GLOBAL LO ADLIB VFLODLIB 

FILEDEF 5 DISK NME DATA 

LOAD ME H ( START 

EXECUTION BEGINS... 

SPECIAL PROJECT : MAXIMUM ENTROPY ALGORITHM 


A MATRIX 2 ROWS 

1 . OOOOOOOD+OO 1. 0000000D+00 
0 . 0000000D+00 1 . OOOOOOOD+OO 


2 COLUMNS 


B MATRIX 
0. OOOOOOOD+OO 
1. OOOOOOOD+OO 


2 ROWS 


1 COLUMNS 


C MATRIX 1 ROWS 

1. OOOOOOOD+OO 0. OOOOOOOD+OO 


2 COLUMNS 


R MATRIX 2 ROWS 

1 . OOOOOOOD+OO 1 . OOOOOOOD+OO 

1 . OOOOOOOD+OO 1. OOOOOOOD+OO 


2 COLUMNS 


R2 MATRIX 
1 . OOOOOOOD+OO 


1 ROWS 


1 COLUMNS 


V MATRIX 2 ROWS 

1. OOOOOOOD+OO 1. OOOOOOOD+OO 
1 . OOOOOOOD+OO 1. OOOOOOOD+OO 


2 COLUMNS 


V2 MATRIX 
1. OOOOOOOD+OO 


1 ROWS 


1 COLUMNS 


B1 MATRIX 2 ROWS 1 COLUMNS 

0. OOOOOOOD+OO 
0. OOOOOOOD+OO 


*** MATRIX F FOR P-RICCATI *** 
8. 0080020D+00 4. 0020000D+00 





*** MATRIX F FOR Q-RICCATI *** 

' 4. 0020000D+00 8 . 0080020D+00 

*** SOLUTION OF LQG P-RICCATI *** 


PROGRAM TO SOLVE CONTINUOUS STEADY-STATE RICCATI EQUATION BY THE NEWTON ALGORITHM 


A MATRIX 2 ROWS 2 COLUMNS 

I . OOOOOOOD+OO 1 . OOOOOOOD+OO 
O.OOOOOOOD+OO 1. OOOOOOOD+OO 

B MATRIX 2 ROWS 1 COLUMNS 

O.OOOOOOOD+OO 
1. OOOOOOOD+OO 


Q MATRIX 2 ROWS 2 COLUMNS 

6. OOOOOOOD+Ol 6 . OOOOOOOD+Ol 
6. OOOOOOOD+Ol 6. OOOOOOOD+Ol 


H IS AN IDENTITY MATRIX 


R MATRIX 1 ROWS 1 COLUMNS 

1. OOOOOOOD+OO 

INITIAL F MATRIX 


F MATRIX 1 ROWS 

8. 0080020D+00 4. 0020000D+00 


2 COLUMNS 


FINAL VALUES OF P AND F AFTER 7 ITERATIONS TO CONVERGE 


P MATRIX 2 ROWS 2 COLUMNS 

2. OOOOOOOD+Ol 1. OOOOOOOD+Ol 
1. OOOOOOOD+Ol 1. OOOOOOOD+Ol 


F MATRIX 1 ROWS 

1. OOOOOOOD+Ol 1. OOOOOOOD+Ol 


2 COLUMNS 


*** SOLUTION OF LQG Q-RICCATI *** 


PROGRAM TO SOLVE CONTINUOUS STEADY-STATE RICCATI EQUATION BY THE NEWTON ALGORITHM 


A MATRIX 

2 ROWS 

2 

COLUMNS 

1. OOOOOOOD+OO 

O.OOOOOOOD+OO 



1. OOOOOOOD+OO 

1. OOOOOOOD+OO 



B MATRIX 

2 ROWS 

1 

COLUMNS 

1. OOOOOOOD+OO 




O.OOOOOOOD+OO 




Q MATRIX 

2 ROWS 

2 

COLUMNS 

6. OOOOOOOD+Ol 

6. OOOOOOOD+Ol 
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6. OOOOOOOD+Ol 6 . 0000000D+01 

H IS AN IDENTITY MATRIX 


R MATRIX 1 ROWS 1 COLUMNS 

1 . OOOOOOOD+OO 

INITIAL F MATRIX 


F MATRIX 1 ROWS 2 COLUMNS 

4. 0020000D+00 8 . 0080020D+00 


FINAL VALUES OF P AND F AFTER 7 ITERATIONS TO CONVERGE 


P MATRIX 2 ROWS 2 COLUMNS 

1. OOOOOOOD+Ol 1. OOOOOOOD+Ol 
1. OOOOOOOD+Ol 2. OOOOOOOD+Ol 

F MATRIX 1 ROWS 2 COLUMNS 

1 . OOOOOOOD+O 1 1 . OOOOOOOD+O 1 

DIF. OF PQ-LYAPUNOV = 1397.87078471772827 

DIF. OF PQ-LYAPUNOV = 0. 568434188608080149E-12 

*** SOLUTION OF ME ALGORITHM *** 

BETA= 0 . OOOOOOOOOOOOOOOOOOE+OO 

*** MATRIX AC *** 

-9. OOOOOOOD+OO 1. OOOOOOOD+OO 

-2. OOOOOOOD+Ol -9. OOOOOOOD+OO 

*** MATRIX F *** 

1. OOOOOOOD+Ol 
1. OOOOOOOD+Ol 

*** MATRIX K *** 

1. OOOOOOOD+Ol 1. OOOOOOOD+Ol 

DIF. OF PQ-LYAPUNOV = 1483.52550453469172 

DIF. OF PQ-LYAPUNOV = 31.1122528653733639 

DIF. OF PQ-LYAPUNOV = 4.03515303082116361 

DIF. OF PQ-LYAPUNOV = 0.141727321507062243 

DIF. OF PQ-LYAPUNOV = 0. 671832436983663683E-01 

DIF. OF PQ-LYAPUNOV = 0. 182016129660951265E-01 

DIF. OF PQ-LYAPUNOV = 0. 233668764548156105E-02 

DIF. OF PQ-LYAPUNOV = 0. 102304770734917838E-03 

*** SOLUTION OF ME ALGORITHM *** 

BETA= 0. 50000000 74505 8059 7E -01 

*** MATRIX AC *** 

-9 . 2759643D+00 1 . OOOOOOOD+OO 

-2. 2199309D+01 -8. 6775519D+00 

*** MATRIX F *** 

1 . 0275964D+01 



1.2521757D+01 


*** MATRIX K *** 

9. 6775519D+00 9. 6775519D+00 


DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

5S 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

S 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

S 

DIF. 

OF 

PQ-LYAPUNOV 

35 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

= 

DIF. 

OF 

PQ-LYAPUNOV 

S 


1549.05227113589928 
84.4184428246941252 
26.0149127163574008 
3.20086419084577756 
2.04724222381747722 
1.86151733248436813 
0.878094194215236712 
0.263785988925747006 
0. 262883900012980121E-01 
0. 268890374742341010E-01 
0 . 2 149898829 15 1197 79E-01 
0. 970914569580827447E-02 
0. 267261225042147998E-02 
0 . 2394 1 32463336973 16E-03 


*** SOLUTION OF ME ALGORITHM *** 
BETA= 0. 999999642372131348E-01 

*** MATRIX AC *** 

-9. 7125436D+00 1 . 0000000D+00 

-2. 4992726D+01 -7 . 325975 1D+00 


*** MATRIX F *** 

1. 0712544D+01 
1. 6666751D+01 

*** MATRIX K *** 

8. 3259751D+00 8. 3259751D+00 


DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV ■ 
DIF. OF PQ-LYAPUNOV » 
DIF. OF PQ-LYAPUNOV - 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV * 
DIF. OF PQ-LYAPUNOV » 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV * 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV « 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 


1665.28167831654781 
159.062705845073140 
71.1380666167275990 
13.3069521641417623 
11.0438398010626315 
16.0436452531334908 
13.0792646555584611 
8.25789695673404367 
4. 14705313050279756 
1.45431039864143941 
0. 465922398768725543E-01 
0.485000639653435428 
0.539609938257910926 
0.402903454531099214 
0.236057926695423248 
0. 106004689502810834 
0. 279536444650148042E-01 
0.9 15 7406695635472 72E-02 
0. 19 72663 140 12836 119E-01 
0. 1829 1252 1599405409E-01 
0. 1204206 79 7 7461 10 17E-01 
0. 652838842739811298E-02 
0. 257274780449279206E-02 
0. 162227934140446450E-03 


*** SOLUTION OF ME ALGORITHM *** 



BETA= 0.149999976158142090 


*** MATRIX AC *** 

- 1 . 0182880D+01 1 . OOOOOOOD+OO 

-2.7716769D+01 -5 . 37 12423D+00 

*** MATRIX F *** 

1. 1182880D+O1 
2. 1345526D+01 

*** MATRIX K *** 

6. 3712423D+00 6. 3712423D+00 


DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV * 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 
DIF. OF PQ-LYAPUNOV = 


2043.32093212088944 
422.622199510942664 
321.264125429695071 
192.278331629577451 
79.8465890746290938 
0.861929903368036321 
45.4660100056273109 
67.0896696260947465 
72.6089619424420221 
68.7181502289284936 
59.9148116939483657 
49.0218397391997769 
37.7723155764265357 
27.2175524990368558 
17.9715890001505159 
10.3474085776692846 
4.43796081232255801 
0.174540652494613369 
2.62192679792053696 
4. 19708560874767045 
4.82445353761380602 
4.77252454138550775 
4.28303268764602763 
3.55564212377225886 
2.74423996702182649 
1.95656871106899644 
1.26023106524792183 
0.689740011255651098 
0.255423312301900296 
0. 497768301354426512E-01 
0.242425045327479438 
0.343552133679565941 
0.376151789676043791 
0.361426173297275000 
0.317537609715657254 
0.258647684084621687 
0. 196205113078065096 
0. 137271073638146390 
0. 858135003701931964E-01 
0. 444365240942943274E-01 
0. 13 7696488600909 106E-01 
0. 759230053478177069E-02 
0. 207745545 176294399E -01 
0. 272035746581309468E-01 
0. 286712760499199248E-01 
0 . 2695695 17 134300440E-01 
0. 232341396629749397E-01 



DIF. OF PQ-LYAPUNOV = 0. 18451 139772594 1606E-01 

' DIF. OF PQ-LYAPUNOV = 0. 139755047 704852586E-01 

DIF. OF PQ-LYAPUNOV = 0. 9728900905 145 1 1030E.-02 

DIF. OF PQ-LYAPUNOV = 0.593533052074235457E-02 

DIF. OF PQ-LYAPUNOV = 0. 2740632 13629460734E-02 

DIF. OF PQ-LYAPUNOV = 0. 648672616364365240E-03 

*** SOLUTION OF ME ALGORITHM *** 

BETA= 0.199999988079071045 

*** MATRIX AC *** 

-1.0741235D+01 1 . OOOOOOOD+OO 

-3. 1813464D+01 -3. 6263966D+00 

*** MATRIX F *** 

1. 1741235D+01 
2. 7187067D+01 

*** MATRIX K *** 

4. 6263966D+00 4 . 6263966D+00 



APPENDIX 2 


Program for M.E. Method 



o n o a a 


C MAIN PROGRAM FOR THE MAXIMUM ENTROPY METHOD 
IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A(10),B(10),C(10),R(10),R1(10),R2(10),V(10), 

& V2(10),B1(10),V1(10), DUMMY( 100) ,FP( 10) , IOP( 3) , 

& AT( 10) ,CT( 10) ,FQ( 10) , H( 10) , P( 10) ,Q( 10) ,PB( 10) , 

& QB( 10) , AS( 10) ,BS( 10) , V2S( 10) ,CS( 10) , BST( 10) , 

& CST( 10) , AST( 10) ,CQ( 10) ,COF( 10),C0P( 10) ,C0P1( 10) , 

& QS( 10) , AQS( 10) ,COQ( 10) ,C0Q1( 10) ,APS( 10), AC 1(10), 

& AC(10),F(10) , AK( 10),UI(10), R2S( 10) , PS( 10) , 

& AP( 10) , AQ( 10) 

DIMENSION NA(2),NB(2),NC(2),NR(2),NR2(2),NV2(2),NB1(2), 

& NV(2),NR1(2) ,NV1(2) ,NCT(2) ,NFP(2) ,NFQ(2) ,NH(2) , 

& NP(2) , NQ(2) ,NAS(2) ,NV2S(2) ,NCS(2) ,NBST(2) , 

& NCST(2),NAST(2),NPS(2),NC0P(2),NAP(2),NAQ(2), 

& NQS(2) ,NAQS(2) ,NC0F(2) ,NCQ(2) ,NAPS(2) ,NAC1(2) , 

& NAC( 2) ,NF(2) ,NK(2) 

LOGICAL IDENT, DISC , FNULL, SYM 

DATA STOL/l.E-4/, ETOL/1 . E-3/ ,EPSA/1 . E-4/.EPSB/1 . E-4/ 

CALL RDTITL 

C INPUT THE MATRICES FOR THE SYSTEM 

CALL READ( 5 , A , NA , B , NB , C , NC , R , NR , R2 , NR2 ) 

CALL READ(3,V,NV,V2,NV2,B1,NB1) 

THETA=60. 

AMU=60. 

CALL SCALE(R,NR,R1 ,NR1 , THETA) 

CALL SCALE(V,NV, V1,NV1 , AMU) 

WRITE(*,*) ' MATRIX Rl' 

CALL PRNT(R1 ,NR1 , 0, 3) 

WRITE(*,*) ' MATRIX Vl' 

CALL PRNT(V1 ,NV1 ,0,3) 

COMPUTE THE F MATRICES FOR P & Q - RICCATI EQUATION 
IOP(1)=0 
I0P(2)=1 
IOP(3)=0 
SCLE=1 . 

CALL CSTAB ( A , NA , B , NB , FP , NFP , I OP , SCLE , DUMMY) 

CALL TRANP(A,NA, AT,NA) 

CALL TRANP(C,NC,CT,NCT) 

CALL CSTAB ( AT, NA , CT, NCT , FQ , NFQ , IOP , SCLE , DUMMY) 

WRITE(*,*) ' *** MATRIX F FOR P-RICCATI ***' 

CALL PRNT(FP,NFP,0,3) 

VIRITE(*,*) ' *** MATRIX F FOR Q-RICCATI ***' 

CALL PRNT(FQ,NFQ,0,3) 

C SOLVE FOR INITIAL P & Q FROM LQG SOLUTION 
I0P(1)=1 
IOP(2)=0 
I0P(3)=0 
IDENT=.TRUE. 

DISC=. FALSE. 

FNULL=. FALSE. 

WRITE(*,*) ' *** SOLUTION OF LQG P-RICCATI ***' 

CALL RICNWT(A,NA,B,NB,H,NH,R1 > NR1,R2,NR2,FP,NFP,P,NP, IOP, 

& I DENT, DISC, FNULL, DUMMY) 

WRITE(*,*) ' *** SOLUTION OF LQG Q-RICCATI ***’ 

CALL RICNWT(AT,NA,CT,NCT,H,NH, V1,NV1 , V2,NV2,FQ,NFQ,Q,NQ, IOP, 
& IDENT, DISC, FNULL, DUMMY) 

C PREPARE THE REQUIRED MATRICES FOR ME ITERATIONS 
CALL NULL(PB.NA) 

CALL NULL(QB,NA) 

CALL EQUATE(A,NA, AS,NAS) 


ME 00010 
ME 00020 
ME 00030 
ME 00040 
ME 00050 
ME 00060 
ME 00070 
ME 00080 
ME 00090 
ME 00100 
ME 00110 
ME 00120 
ME 00130 
ME 00140 
ME 00150 
ME 00160 
ME 00170 
ME 00180 
ME 00190 
ME 00200 
ME 00210 
ME 00220 
ME 00230 
ME 00240 
ME 00250 
ME 00260 
ME 00270 
ME 00280 
ME 00290 
ME 00300 
ME 00310 
ME 00320 
ME 00330 
ME 00340 
ME 00350 
ME 00360 
ME 00370 
ME 00380 
ME 00390 
ME 00400 
ME 00410 
ME 00420 
ME 00430 
ME 00440 
ME 00450 
ME 00460 
ME 00470 
ME 00480 
ME 00490 
ME 00500 
ME 00510 
ME 00520 
ME 00530 
ME 00540 
ME 00550 
ME 00560 
ME 00570 
ME 00580 
ME 00590 
ME 00600 
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_ 

CALL EQUATE (B.NB.BS.NBS) 

ME 

00610 


CALL EQUATE (V2.NV2 ,V2S,NV2S) 

ME 

00620 


CALL EQUATE ( C , NC , CS , NCS ) 

ME 

00630 


CALL TRANP(BS ,NBS, BST, NBST) 

ME 

00640 


CALL TRANP(CS , NCS , CST, NCST) 

ME 

00650 


CALL TRANP( AS , NAS , AST, NAST) 

ME 

00660 


CALL UNITY(UI.NA) 

ME 

00670 


S=-l. 

ME 

00680 


DO 300 IK=1 ,5 

ME 

00690 


BETA=. 05*( IK- 1 ) 

ME 

00700 


B1(2)=BETA 

ME 

00710 

c 

BEGIN ITERATIONS 

ME 

00720 


PQTEMP=0 . 

ME 

00730 

5 

K=1 

ME 

00740 


PTNORM=0 . 

ME 

00750 


QTN0RM=0. 

ME 

00760 


PLTN0R=0 . 

ME 

00770 


QLTN0R=0 . 

ME 

00780 

C 

COMPUTE COEFFICIENTS FOR P-RICCATI 

ME 

00790 


1=1 

ME 

00800 

10 CALL SUB 1 2 ( R2 , NR2 , B 1 , NB 1 , P , NP , PB , NA , R2S , NR2 ) 

ME 

00810 

C 

SOLVE P-RICCATI 

ME 

00820 


IOP(1)=0 

ME 

00830 


IOP(2)=0 

ME 

00840 


I0P(3)=0 

ME 

00850 


IDENT=. TRUE . 

ME 

00860 


DISC=.FALSE. 

ME 

00870 


FNULL=. FALSE. 

ME 

00880 

C 

WRITE(*,*) ’ *** SOLUTION OF P-RICCATI ***' 

ME 

00890 


CALL RICNWT(AS,NA,BS,NBS,H,NH,R1,NR1,R2S,NR2,FP,NFP,P,NP, IOP, 

ME 

00900 


& IDENT, DISC, FNULL, DUMMY) 

ME 

00910 

C 

TEST FOR CONVERGENCE OF P - RICCATI SOLUTION 

ME 

00920 


I0PT=2 

ME 

00930 


M1=NP( 1) 

ME 

00940 


CALL N0RMS(M1,M1,M1,P, I0PT,PN0RM) 

ME 

00950 


DIF=DABS(PNORM-PTNORM) 

ME 

00960 

C 

WRITE(*,*) 1 DIF. OF P-RICCATI = ' , DIF 

ME 

00970 


IF(DIF.LE.STOL) THEN 

ME 

00980 


GO TO 20 

ME 

00990 


ELSE 

ME 

01000 


PTNORM=PNORM 

ME 

01010 


1=1+1 

ME 

01020 


IF(I.GE.SOO) GO TO 200 

ME 

01030 


GO TO 10 

ME 

01040 


END IF 

ME 

01050 

C 

COMPUTES COEFFICIENT FOR Q - RICCATI EQUATION 

ME 

01060 

20 

i J=1 

ME 

01070 

25 

CALL SUB12(R2,NR2,B1,NB1,P,NP,PB,NA,R2S,NR2) 

ME 

01080 


CALL MULT ( BST , NBST , P , NP , PS , NPS ) 

ME 

01090 


CALL SUB 13(B1,NB1,R2S,NR2,PS, NPS , QB , NA , CQ , NCQ) 

ME 

01100 


CALL ADD( V 1 , NV 1 , CQ , NCQ , COF , NCOF ) 

ME 

onio 

C 

SOLVE FOR Q-RICCATI 

ME 

01120 

C 

VRITE(*,+) ' *** SOLUTION OF Q-RICCATI EQ. ' 

ME 

01130 


CALL R I CNWT( AST , NAST , CST , NCST , H , NH , COF , NCOF , V2S , NV2S , FQ , NFQ , 

ME 

01140 


& Q,NQ, IOP, IDENT, DISC, FNULL, DUMMY) 

ME 

01150 

C 

TEST FOR CONVERGENCE OF Q-RICCATI 

ME 

01160 


N1=NQ( 1). 

ME 

01170 


CALL N0RMS(N1,N1,N1,Q,I0PT,QN0RM) 

ME 

01180 


D I F=DAB S ( QNORM - QTNORM ) 

ME 

01190 

C 

WRITE(*,*) 1 DIF. OF Q-RICCATI = ' , DIF 

ME 

01200 
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IF(DIF.LE.STOL)THEN 
GO TO 30 

ELSE 

QTNORM=QNORM 

J=J+1 

IF( J.GE. 500) GO TO 200 
GO TO 25 
END IF 

C COMPUTE COEFFICIENTS FOR P-LYAPUNOV EQUATION 
30 11=1 

35 CALL SUB 1 2 ( R2 , NR2 , B 1 , NB 1 , P , NP , PB , NA , R2S , NR2 ) 

ITYPE=1 

CALL SUB5C I TYPE , UI , NA, P , NP, BS , NBS , R2S ,NR2 , COP, NCOP) 

C WRITE(*,*) ' *** MATRIX C OF P-LYAPUNOV ***' 

C CALL PRNT(COP,NCOP, 0, 3) 

CALL SCALE (COP, NCOP, COP 1, NCOP, S) 

CALL MULT(Q,NQ,CST,NCST,QS,NQS) 

CALL SUB8 ( AS , NAS , QS , NQS , V2S , NV2S ,CS , NCS , AQS , NAQS) 

C SOLVE P-LYAPUNOV EQUATION 
I0PL=0 
SYM=. TRUE. 

C WRITE(*,*) ' *** SOLUTION P-LYAPUNOV EQ. ***' 

CALL BARSTW( AQS , NAQS , AQ , NAQ , COP 1 , NCOP , IOPL, SYM, EPSA , EPSB , DUMMY) 
CALL EQUATE ( COP 1, NCOP, PB.NA) 

C TEST FOR CONVERGENCE OF P-LYAPUNOV 

CALL NORMS ( Ml, Ml, Ml, PB, IOPT.PLNORM) 

DIF=DABS(PLTNOR-PLNORM) 

C WRITE(*,*) ' DIF. OF P-LYAPUNOV =',DIF 
IF(DIF.LE.STOL) THEN 
GO TO 40 
ELSE 

PLTNOR=PLNORM 
IF(Il.GE.SOO) GO TO 200 
GO TO 35 
END IF 

C COMPUTE COEFFICIENTS FOR Q-LYAPUNOV EQUATION 


C40 

Jl=l 

C45 

ITYPE=2 

40 

ITYPE=2 


CALL SUBS ( I TYPE , UI , NA , Q , NQ , CS , NCS , V2S , NV2S , COQ, NCOQ) 

C WRITE(*,*) ' *** MATRIX C OF Q-LYAPUNOV ***’ 

C CALL PRNT(COQ, NCOQ, 0,3) 

CALL SCALE(COQ,NCOQ,COQl,NCOQ,S) 

CALL SUB 1 2 ( R2 , NR2 , B 1 , NB 1 , P , NP , PB , NA , R2S , NR2 ) 

CALL MULT( BST , NBST, P , NP , PS , NPS ) 

CALL SUB8 ( AS , NAS , BS , NBS , R2S , NR2 , PS , NPS , APS , NAPS ) 

C SOLVE Q-LYAPUNOV EQUATION 

C WRITE(*,*) ' *** SOLUTION OF Q-LYAPUNOV ***' 

CALL BARSTVC APS , NAPS , AP , NAP , COQ 1 , NCOQ , IOPL, SYM, EPSA , EPSB , DUMMY ) 
CALL EQUATE ( COQ 1, NCOQ, QB.NA) 

C TEST FOR CONVERGENCE OF Q-LYAPUNOV 

CALL NORMS (Ml, Ml, Ml, QB, IOPT.QLNORM) 

C DIF=DABS(QLNORM-QLTNOR) 

C WRITE(*,*) ' DIF. OF Q-LYAPUNOV =’,DIF 

C IF(DIF.LE.STOL) THEN 

C GO TO 50 

C ELSE 

C QLTNOR=QLNORM 

C J1=J1+1 

C IF( J1 . GE. 500) GO TO 200 


ME 01210 
ME 01220 
ME 01230 
ME 01240 
ME 01250 
ME 01260 
ME 01270 
ME 01280 
ME 01290 
ME 01300 
ME 01310 
ME 01320 
ME 01330 
ME 01340 
ME 01350 
ME 01360 
ME 01370 
ME 01380 
ME 01390 
ME 01400 
ME 01410 
ME 01420 
ME 01430 
ME 01440 
ME 01450 
ME 01460 
ME 01470 
ME 01480 
ME 01490 
ME 01500 
ME 01510 
ME 01520 
ME 01530 
ME 01540 
ME 01550 
ME 01560 
ME 01570 
ME 01580 
ME 01590 
ME 01600 
ME 01610 
ME 01620 
ME 01630 
ME 01640 
ME 01650 
ME 01660 
ME 01670 
ME 01680 
ME 01690 
ME 01700 
ME 01710 
ME 01720 
ME 01730 
ME 01740 
ME 01750 
ME 01760 
ME 01770 
ME 01780 
ME 01790 
ME 01800 



C GO TO 45 

C END IF 

C TEST FOR CONVERGENCE OF. ME SOLUTION 
50 PQNORM=PLNORM+QLNORM 

DIF=DABS ( PQTEMP - PQNORM ) 

WRITE(*,*) ' DIF. OF PQ- LYAPUNOV =', DIF 
IF(DIF. LE.ETOL) THEN 
GO TO 60 
ELSE 

PQTEMP=PQNORM 
IF(K. GE . 50) GO TO 200 
GO TO 10 
END IF 

C COMPUTE COMPENSATER MATRICES 
C COMPUTE AC 

60 CALL SUBS ( AS , NAS , QS , NQS , V2 S , NV2 S , CS , NCS , AC 1 , NAC 1 ) 
CALL SUB 1 2 ( R2 , NR2 , B 1 , NB 1 , P , NP , PB , NA , R2S , NR2 ) 

CALL SUB 8 ( AC 1 , NAC 1 , B S , NB S , R2 S , NR2 , PS , NPS , AC , NAC ) 
WRITE(*,*) ' ' 

WRITE(*,*) ' *** SOLUTION OF ME ALGORITHM ***' 
WRITE(*,*) ' BETA=', BETA 
WRITE(*,*) ' *** MATRIX AC ***' 

CALL PRNT( AC , NAC ,0,3) 

C COMPUTE F 

ITYPE=2 

CALL SUB9 ( I TYPE , Q , NQ , CS , NCS , V2S , NV2S , F , NF) 
VRITE(*,*) ' *** MATRIX F ***' 

CALL PRNT(F,NF,0,3) 

C COMPUTE K 

ITYPE=1 

CALL SUB9(ITYPE,R2S,NR2,BS,NBS,P,NP > AK,NK) 
WRITE(*,*) ' *** MATRIX K ***' 

CALL PRNT(AK,NK, 0, 3) 

300 CONTINUE 
200 STOP 
END 

C ****** SUBROUTINE SUB1 

SUBROUTINE SUB1(B,NB,C,NC,D,ND,A,NA) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A(50) , B(50) ,C(50) ,D(50) , BC(50) 

DIMENSION NA(2),NB(2),NC(2),ND(2),NBC(2) 

CALL MULT( B , NB , C , NC , BC , NBC ) 

CALL MULT(BC,NBC,D,ND,A,NA) 

RETURN 

END 

C ****** SUBROUTINE SUB5 

SUBROUTINE SUBSCITYPE.B.NB.C.NC.D.ND.E.NE.A.NA) 
IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A(50),B(50),C(50),D(50),E(50), 

& DT(50),F(50) ,FT(50) ,EI(50) ,BT(50) 

DIMENSION NA(2),NB(2),NC(2),ND(2),NE(2),NBT(2), 

& NDT(2),NF(2),NFT(2) 

CALL TRANP(B,NB, BT,NBT) 

IF(ITYPE.EQ. 1) CALL SUB1(BT,NBT,C,NC,D,ND,F,NF) 
IF(ITYPE.EQ. 2) CALL SUBKD.ND.C.NC.B.NB.F.NF) 

CALL TRANP(F,NF,FT,NFT) 

CALL UNITY(EI.NE) 

N=NE( 1) 

NR=NE(2) 

CALL GAUSEL(N,N,E,NR,EI, IERR) 
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IF( ITYPE. EQ. 1) CALL SUB 1(F, NF , El ,NE , FT.NFT, A,NA) 

IF( ITYPE.EQ. 2) CALL SUB1(FT,NFT,EI ,NE , F,NF, A,NA) 

RETURN 

END 

C ***** SUBROUTINE SUB8 

SUBROUTINE SUB8 ( B , NB , C , NC , D , ND , E , NE , A , NA) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION B(50),C(50),D(50),E(50),A(50),F(50) 

DIMENSION NB(2),NC(2),ND(2) ,NE(2) , NA(2),NF(2) 

CALL UNITY(DI.ND) 

N=ND( 1) 

NR=ND(2) 

CALL GAUSEL(N,N,D,NR,DI, IERR) 

CALL SUB 1 (C , NC , DI , ND , E , NE , F , NF) 

CALL SUBT( B , NB , F , NF , A , NA) 

RETURN 

END 

C ****** SUBROUTINE SUB 9 

SUBROUTINE SUB9( ITYPE.B.NB.C.NC.D.ND.A.NA) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A(50) , B(50) , C(50) ,D(50) , BI(50) ,CI(50) ,DI(50) ,CT(50) 
DIMENSION NA(2) , NB(2) ,NC(2) ,ND(2) ,NCT(2) 

IF( ITYPE. EQ. 1) THEN 
CALL UNITY(BI.NB) 

N=NB( 1) 

NR=NB ( 2 ) 

CALL GAUSEL(N,N,B,NR,BI, IERR) 

ELSE 

CALL UNITY(DI.ND) 

N=ND( 1) 

NR=ND(2) 

CALL GAUSEL(N,N,D,NR,DI, IERR) 

END IF 

CALL TRANP(C,NC,CT,NCT) 

IF( ITYPE. EQ. 1) CALL SUB1(BI,NB,CT,NCT,D,ND,A,NA) 

IF( ITYPE.EQ. 2) CALL SUB1(B,NB,CT,NCT,DI ,ND,A,NA) 

RETURN 

END 

C ***** SUBROUTINE SUB 12 

SUBROUTINE SUB12(A,NA,B,NB,C,NC,D,ND,E,NE) 

IMPLICIT REAL* 8 (A-H.O-Z) 

DIMENSION A(50) , B(50) ,C(50) ,D(50) ,E(50) ,BT(50) ,CD(50) ,TEMP(50) 
DIMENSION NA(2) , NB(2) ,NC(2) ,ND(2) ,NE(2) ,NBT(2) ,NCD(2) 

CALL TRANP(B,NB,BT,NBT) 

CALL ADD(C,NC,D,ND,CD,NCD) 

CALL SUB 1 ( BT , NBT , CD , NCD , B , NB , TEMP , NA ) 

CALL ADD(A,NA,TEMP,NA,E,NE) 

RETURN 

END 

C**** SUBROUTINE SUB 13 

SUBROUTINE SUB13(A,NA,B,NB,C,NC,D,ND,E,NE) 

DIMENSION A(50),B(50),C(50),D(50),E(50),BI(50) ,TEMP(50) ,TT(50) 
DIMENSION NA(2) , NB(2) , NC(2) ,ND( 2) ,NE(2) ,NT(2) ,NTT(2) 

CALL UNITY(BI.NB) 

N=NB( 1) 

NR=NB ( 2 ) 

CALL GAUSEL(N,N,B,NR,BI, IERR) 

CALL SUB 1 ( A , NA , B I , NB , C , NC , TEMP , NT ) 

CALL TRANP(TEMP, NT.TT.NTT) 

CALL SUBKTEMP.NT.D.ND.TT.NTT.E.NE) 
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RETURN 

END 
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C MAIN PROGRAM FOR THE OPTIMAL PROJECTION METHOD 
IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION A(49),B(14),C(21),R1(49),R2(4),V1(49), 

& V2(9),P(49),Q(49),UI(49),TAUO(49),Cl(49), 

& IOP(3),F(49),C3(49),CT(21),C5(49),C6(49), 

& C12(21) ,AQC(49) ,AQ(49) ,AQT(49) ,BX(49) , 

& C8(49) ,C9(49) ,C13( 14) , AP(49) ,APC(49) ,ER(50) , 

& EI(57),V(49) ,TAU(49) ,C1 1(49) ,C14(49) ,GA(49) , 

& G(49) ,GT(49) ,AC(49) ,FC(49) ,RKC(49) ,H(49) , 

& FP(49) ,FQ(49) ,DUMMY(500) ,AT(49) ,APT(49) , 

& VTK( 98) ,QP(49) ,CK( 14) ,CF(21) ,CFC(49) ,R2N(49) , 

& BCK(49) ,ACF(49) ,CA(49) ,AP1(49) ,AQ1(49) ,VI(49) 

DIMENSION NA(2),NB(2),NC(2),NR2(2),NV2(2),NF(2),NCT(2), 

& NBX(2) ,NC13(2) ,NGA(2) ,NG(2) ,NGT(2) ,NAC(2) , 

& NRKC(2),NH(2),NFP(2),NFQ(2),NR1(2),NV1(2), 

& NP(2) ,NQ(2) ,NCK(2) ,NCF(2) ,NAP1(2) ,NAQ1(2) ,NC12(2) 

LOGICAL I DENT, DISC, FNULL, SYM 

DATA STOL/l.E-4/, ETOL/l.E-3/,EPSA/l.E-6/,EPSB/l.E-6/ 

CALL RDTITL 

C WRITE(*,*) ' INPUT THE ORDER TO BE REDUCED ' 

C READ(*,*) NCR 

NCR=4 

C INPUT THE MATRICES FOR THE SYSTEM 

CALL READ ( 5 , A , NA , B , NB , C , NC , R 1 , NR 1 , R 2 , NR2 ) 

CALL READ(2,V1,NV1,V2,NV2) 

C R2(2)=l . 

C R2(3)=2. 

WRITE(6,*) ' *** NORMAL R2' 

CALL NORMAL(R2,NR2,R2N,NR2) 

CALL PRNT(R2N,NR2,0,3) 

C COMPOTE THE F MATRICES FOR P & Q - RICCATI EQUATION 
IOP(1)=0 
IOP(2)=l 
IOP(3)=0 
SCLE=1 

CALL CSTAB(A,NA,B,NB,FP,NFP,IOP,SCLE, DUMMY) 

WRITE(6,*) ' MATRIX F' 

CALL TRANP(A,NA,AT,NA) 

CALL TRANP(C,NC,CT,NCT) 

CALL CSTAB(AT,NA,CT,NCT,FQ,NFQ,IOP,SCLE, DUMMY) 

C SOLVE FOR INITIAL P & Q FROM LQG SOLUTION 
IOP(1)=0 
IOP(2)=0 
IOP(3)=0 
IDENT= . TRUE . 

DISC=. FALSE. 

FNULL=. FALSE. 

WRITE(6,*) ' RICCATI' 

CALL RICNWT(A,NA,B,NB,H,NH,R1,NR1,R2,NR2,FP,NFP,P,NP,I0P, 

& IDENT, DISC, FNULL, DUMMY) 

WSITE(6,*) ' Q RICCATI’ 

CALL RICNWT(AT,NA ,CT,NCT,H,NH,V1,NV1 ,V2,NV2,FQ,NFQ,Q,NQ , IOP , 
& I DENT, DISC, FNULL, DUMMY) 

C COMPUTE THE COMPENSATOR MATRICES FOR LQG 
CALL SUB9( 1 ,R2 ,NR2 ,B ,NB ,P,NP,CK,NCK) 

CALL SUB9 ( 2 , Q , NQ , C , NC , V2 , NV2 , CF , NCF ) 

CALL MULT(CF,NCF,C,NC,CFC,NA) 

CALL MULT(B,NB,CK,NCK,BCK,NA) 

CALL SUBT ( A , NA , CFC , NA , ACF , NA ) 
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CALL SUBT(ACF,NA , BCK ,NA ,CA ,NA) 

WRITE(6,*) ' K MATRIX FOR COMPENSATOR’ 

CALL PRNT(CK,NCK,0 ,3) 

WRITE(6 ,*) ' F MATRIX FOR COMPENSATOR' 

CALL PRNT(CF,NCF,0,3) 

WRITE(6,*) ' AC MATRIX FOR COMPENSATOR' 

CALL PRNT(CA,NA,0,3) 

C COMPUTES MATRIX NORM P & Q SOLUTIONS 
M1=NA( 1) 

N1=NA(2) 

IOPT=2 

WRITE(6,*) ' NOW A' 

C CALL NORMS(Ml,Ml,Nl,P,IOPT,PNORM) 

WRITE(6,*) ' NOW B' 

C CALL N0RMS(M1,M1,N1,Q, IOPT,QNORM) 

CALL UNITY(UI.NA) 

CALL NULL(TAUO.NA) 

C BEGIN ITERATIONS FOR OPTIMAL PROJECTION ALGORITHM 
K=1 

5 1=1 

PNORM=0 . 

C COMPUTES COEFFICIENT FOR P - RICCATI EQUATION 
10 ITYPE=1 

WRITE(6,*) ' NOW C' 

CALL SUB5 ( ITYPE , TAUO , NA , P , NA , B , NB , R2 , NR2 , C 1 , NA ) 
WRITE(6,*) ' NOW D' 

CALL ADD(R1,NR1,C1,NA,C1,NA) 

WRITE(*,*) ' NOW E' 

C SOLVES FOR P - RICCATI EQUATION 
IOP(1)*0 
IOP(2)=0 
IOP(3)*0 
I DENT*. TRUE. 

DISC*. FALSE. 

FNULL*. FALSE. 

CALL RICNWT(A,NA,B,NB,H,NH,C1,NA,R2,NR2,FP,NFP,P,NP,I0P, 
& I DENT, DISC, FNULL, DUMMY) 

WRITE(*,*) ' PASS P-RICCATI * 

C TEST FOR CONVERGENCE OF P - RICCATI SOLUTION 
I0PT=2 

CALL N0RMS(M1,M1,N1,P,I0PT,PTN0RM) 

DIF=DABS( PNORM - PTNORM ) 

WRITE(*,*) ' DIF*', DIF 
IF(DIF.LE.STOL) THEN 
GO TO 20 
ELSE 

PNORM=PTNORM 

1=1+1 

IF(I.GE.IOOO) GO TO 200 
GO TO 10 
END IF 
20 J=1 

QN0RM=0 . 

C COMPUTES COEFFICIENT FOR Q - RICCATI EQUATION 
WRITE(*,*) ’ NOW ONE' 

30 ITYPE=2 

CALL SUB5( ITYPE, TAUO, NA,Q,NA,C,NC,V2,NV2,C3,NA) 

CALL ADD(V1,NA,C3,NA,C3,NA) 

C SOLVES FOR Q - RICCATI EQUATION 
WRITE(*,*) ' NOW Q' 
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CALL RICNWT(AT,NA ,CT ,NCT,H ,NH ,C3 ,NA , V2 ,NV2 ,FQ ,NFQ ,Q,NQ , IOP , 
& IDENT .DISC, FNULL , DUHMY ) 

C TEST FOR CONVERGENCE OF Q - RICCATI SOLUTION 
WRITE(*,*) 1 NORMS' 

CALL NORMS (M1,M1,N1,Q, IOPT , QTNORM) 

DIF=DABS(QNORM-QTNORM) 

WRITE(*,*) ' DIFQ=',DIF 
IF(DIF.LE.STOL) THEN 
GO TO 40 
ELSE 

QNORM=QTNORM 

J=J+1 

IF(J.GE.IOOO) GO TO 200 
WRITE(*,*) ' GO TO 30’ 

GO TO 30 
END IF 

C COMPUTE COEFFICIENTS FOR P-LYAPUNOV EQUATION 
40 ITYPE=1 

WRITE(*,*) ’ NOW TWO’ 

CALL SUB5(ITYPE,UI,NA,P,NA,B,NB,R2,NR2,C5,NA) 

WRITE(*,*) ’ NOW 3’ 

CALL SUBS (ITYPE , TAUO , NA , P , NA , B , NB , R2 , NR2 , C6 , NA ) 

WRITE(*,*) ’ NOW 4’ 

CALL SUBT(C6,NA,C5 ,NA,C6,NA) 

ITYPE=2 

CALL SUB9( ITYPE, Q,NA,C,NC,V2,NV2,C12,NC12) 

WRITE(*,*) 1 N0W5 ' 

CALL MULT(C12,NC12,C,NC,AQC,NA) 

WRITE(*,*) * N0W6 1 

CALL SUBT(A,NA,AQC,NA,AQ,NA) 

WRITE(*,*) ’ AQ BARSTW - P* 

CALL PRNT(AQ,NA,0,3) 

C SOLVE FOR P - LYAPUNOV EQUATION 
IOPL=l 
SYM= . TRUE . 

CALL TRANP(AQ,NA,AQT,NA) 

WRITE(*,*) ’ NOW7 1 

CALL BARSTW(AQT,NA,AQ1,NAQ1,C6,NA,I0PL,SYM,EPSA,EPSB, DUMMY) 
C COMPUTE COEFFICIENTS FOR Q - RICCATI EQUATION 
ITYPE=2 

WRITE(*,*) ’Ql’ 

CALL SUB5( ITYPE, UI,NA,Q,NA,C,NC,V2,NV2,C8,NA) 

WRITE(*,*) ’ Q2' 

CALL SUB5( ITYPE, TAUO, NA,Q,NA,C,NC,V2,NV2,C9,NA) 

WRITE(*,*) ' Q3' 

CALL SUBT(C9,NA,C8,NA,C9,NA) 

ITYPE=1 

CALL SUB9 ( ITYPE, R2,NR2,B,NB,P,NA,C13,NC 13) 

CALL MULT(B,NB,C13,NC13,APC,NA) 

WRITE(*,*) ' Q4' 

CALL SUBT(A,NA,APC,NA,AP,NA) 

WRITE(*,*) ' AP BARSTW - Q' 

CALL PRNT(AP,NA,0,3) 

C SOLVES FOR Q - LYAPUNOV EQUATION 
WRITE(*,*) ' WRITE' 

PATT TPANPfAP NA APT NA'l 

CALL B ARSTW ( AP , NA , AP 1 , NAP 1 , C9 , NA , IOPL , SYM , EPSA , EPSB , DUMMY ) 

C TEST FOR CONVERGENCE OF P & Q - LYAPUNOV SOLUTIONS 
CALL MULT(C9,NA,C6,NA,QP,NA) 

WRITE(*,*) ' *** MATRIX QP ***’ 
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CALL PRNTCQP ,NA ,0,3) 

* C COMPUTE EIGENVALUES AND EIGENVECTORS OF MATRIX QP 
N=NA(1) 

ISV=N 

ILV=0 

CALL EIGEN(N,N,QP,ER,EI,ISV,ILV,V,WK,IERR) 

WRITE(*,*) ' ISV =’,ISV 
WRITE(*,*) ' ILV ■' , ILV 
WRITE(*,*) ' IERR =',IERR 

C CHECK IF EIGENVALUES ARE ARRANGED IN INCREASING OR DECREASIG ORDER 
CALL LNCNT(4) 

PRINT 650 

650 FORMATC//,' EIGENVALUES OF QP’,//) 

675 F0RMAT( 10X.2D16 .8) 

CALL LNCNT(N) 

DO 700 11=1, N 

PRINT 675,ER(I1),EI(I1) 

700 CONTINUE 

WRITE(*,*) ' EIGENVECTOR OF QP WITH NAMDA INCREASING ORDER’ 
CALL PRNT(V,NA,0,3) 

N=NA(1) 

NU=N-NCR 

ND=NU+1 

RA=ER(NU)/ER(ND) 

RATIO=DABS(RA) 

WRITE(*,*) ' RATIO=\ RATIO 
IF ( RAT I 0 . LT . ETOL ) THEN 
GO TO 50 
ELSE 
K=K+1 

IF(K.GE.500) GO TO 200 
C FORM NEW TAU 
C CALL UNITY(VI.NA) 

C N=NA(1) 

C NR=NA(2) 

C CALL GAUSEL(N,N,V,NR,VI,IERR) 

C CALL FOMTAU ( V , NA , NCR , TAU , NA ) 

CALL CONTAU (NCR,VI,NA,TAU,NA) 

CALL SUBT(UI ,NA,TAU,NA,TAUO,NA) 

WRITEC*,*) ' TAU' 

CALL PRNT(TAU,NA,0,3) 

WRITE(*,*) ' TAUO' 

CALL PRNT( TAUO, NA, 0,3) 

WRITEC*,*) ' GO TO 5' 

GO TO 5 
END IF 

50 CALL SUBT(AQ,NA,APC,NA,C11,NA) 

CALL SUB1(C12,NC,D,ND,C13,NC13,C14,NA) 

CALL ADD(C11,NA,C14,NA,C14,NA) 

C FORM GAMMA AND G 

C 

C 

C 

C 

C CALL SUB 1 (GA , NGA , C 14 , NA ,GT , NG , AC , NAC ) 

C PRINT AC 

C CALL MULT(GA,NGA,C12,NC,FC,NFC) 

C PRINT FC 

C CALL MULT(C13,NC13,GT,NG,RKC,NRKC) 

C PRINT KC 
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OP 02300 
OP 02310 
OP 02320 
OP 02330 
OP 02340 
OP 02350 
OP 02360 
OP 02370 
OP 02380 
OP 02390 
OP 02400 
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200 STOP 
END 

C ****** SUBROUTINE SUB1 

SUBROUTINE SUB1(B,NB,C,NC,D,ND,A,NA) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A(*) ,B(*) ,C(*) ,D(*) ,BC(49) 

DIMENSION NA(2) ,NB(2) ,NC(2) ,ND(2) ,NBC(2) 

CALL MULT ( B , NB , C , NC , BC , NBC ) 

CALL MULT(BC ,tyBC,D,ND,A,NA) 

RETURN 

END 

C ****** SUBROUTINE SUBS 

SUBROUTINE SUB5 ( I TYPE , B , NB , C , NC , D , ND , E , NE , A , NA) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A(50) ,B(*) ,C(*) ,D(*) ,E(*) , 

& DT(50),F(50),FT(50),EI(50),BT(50) 

DIMENSION NA(2),NB(2),NC(2),ND(2),NE(2),NBT(2), 

& NDT(2),NF(2),NFT(2) 

CALL TRANP(B,NB,BT,NBT) 

IF(ITYPE.EQ.l) CALL SUB1(BT,NBT,C,NC,D,ND,F,NF) 

IF(ITYPE.EQ.2) CALL SUB1(D,ND,C,NC,BT,NBT,F,NF) 

CALL TRANP(F,NF,FT,NFT) 

CALL UNITY(EI.NE) 

N=NE( 1) 

NR=NE(2) 

CALL GAUSEL(N,N,E,NR,EI,IERR) 

IF(ITYPE.EQ.l) CALL SUB1(F,NF,EI,NE,FT,NFT,A,NA) 
IF(ITYPE.EQ.2) CALL SUB1(FT,NFT,EI ,NE,F,NF,A,NA) 

RETURN 

END 

C ****** SOUROUTINE SUB9 

SUBROUTINE SUB9(ITYPE,B,NB,C,NC,D,ND,A,NA) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A(50),B(50),C(50),D(50),BI(50),CI(50),DI(50),CT(50) 
DIMENSION NA(2) ,NB(2) ,NC(2) ,ND(2) ,NCT(2) 

IF(ITYPE.EQ.l) THEN 
CALL UNITY(BI.NB) 

N=NB(1) 

NR=NB(2) 

CALL GAUSEL(N,N,B > NR,BI, IERR) 

ELSE 

CALL UNITY(DI.ND) 

N=ND( 1) 

NR=ND(2) 

CALL GAUSEL(N,N,D,NR,DI,IERR) 

END IF 

CALL TRANP(C,NC,CT,NCT) 

IF(ITYPE.EQ.l) CALL SUB1(BI,NB,CT,NCT,D,ND,A,NA) 
IF(ITYPE.EQ.2) CALL SUB1(B,NB,CT,NCT,DI,ND,A,NA) 

RETURN 

END 

C ***** SUBROUTINE FOMTAU 

SUBROUT I NE FOMT AU ( V , NV , NCR , TAU , NA ) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION V(50) ,TAU(50) ,VI(50) ,SUM(50) ,VKT(50) ,UK(50) ,UV(50) 
DIMENSION NV(2) ,NA(2) ,NVKT(2) ,NUK(2) 

CALL UNITY(VI.NV) 

N=NV( 1) 

NR=NV( 2) 

CALL GAUSEL(N,N,V,NR,VI , IERR) 


OP 02410 
OP 02420 
OP 02430 
OP 02440 
OP 02450 
OP 02460 
OP 02470 
OP 02480 
OP 02490 
OP 02500 
OP 02510 
OP 02520 
OP 02530 
OP 02540 
OP 02550 
OP 02560 
OP 02570 
OP 02580 
OP 02590 
OP 02600 
OP 02610 
OP 02620 
OP 02630 
OP 02640 
OP 02650 
OP 02660 
OP 02670 
OP 02680 
OP 02690 
OP 02700 
OP 02710 
OP 02720 
OP 02730 
OP 02740 
OP 02750 
OP 02760 
OP 02770 
OP 02780 
OP 02790 
OP 02800 
OP 02810 
OP 02820 
OP 02830 
OP 02840 
OP 02850 
OP 02860 
OP 02870 
OP 02880 
OP 02890 
OP 02900 
OP 02910 
OP 02920 
OP 02930 
OP 02940 
OP 02950 
OP 02960 
OP 02970 
OP 02980 
OP 02990 
OP 03000 
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CALL NULL ( SUM, NV) 

DO 10 K=1 ,NCR 

CALL UKVKT ( K , V , NV , V I , VKT , NVKT , UK , NUK ) 

CALL MULT (UK, NUK, VKT, NVKT, UV,NV) 

CALL ADD( SUM , NV , UV , NV , SUM , NV ) 

10 CONTINUE 

CALL EQUATE ( SUM, NV,TAU,NA) 

RETURN 

END 

C ****** SUBROUTINE UKVKT 

SUBROUTINE UKVKT ( K , V , NV , V I , VKT , NVKT , UK , NUK ) 

IMPLICIT REAL* 8 (A-H.O-Z) 

DIMENSION V(50) ,VI(50) ,VKT(50) ,UK(50) 

DIMENSION NV(2),NVKT(2),NUK(2) 

N=NV( 1) 

L=1+(K-1)*N 
DO 10 1=1, N 
JV=K+(I-1)*N 
VKT( I)=V( JV) 

JU=L+(I-1) 

UK(I)=VI(JU) 

10 CONTINUE 
NVKT ( 1 ) = 1 
NVKT(2)=N 
NUK ( 1 ) =N 
NUK(2)=1 
RETURN 
END 

C ****** SUBROUTINE CONTAU 

SUBROUTINE CONTAU(NCR , PI ,NA , TAU , NTAU) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION PI (49) ,TAU(49) ,PSI(49) ,EI(49),NA(2) ,NTAU(2) ,PN(49) 
C CONSTRUCT PSI FROM PI 

CALL PSICON(PI,NA,PSI,NA) 

WRITE(*,*) ' EIGENVECTOR OF QP WITH NAMDA DECREASING ORDER' 
CALL PRNT(PSI,NA,0,3) 

C CONSTRUCT MATRIX (INC,0) 

CALL NORMAL(PSI,NA,PN,NA) 

WRITE(*,*) ' NORMALIZED EIGENVECTOR' 

CALL PRNT(PN,NA,0,3) 

CALL NULL(EI.NA) 

N=NA(1) 

N1=N+1 

DO 10 1=1, NCR 
K=1+(I-1)*N1 
EI(K)=1 
10 CONTINUE 

WRITE(*,*) • MATRIX (INC, 0)' 

CALL PRNT(EI ,NA,0 ,3) 

C COMPUTES TAU 
ITYPE=2 

CALL SUB9(ITYPE,PN,NA,EI,NA,PN,NA,TAU,NA) 

RETURN 

END 

C ****** SUBROUTINE PS ICON 

SUBROUTINE PSICON(PI,NA,PSI,NPSI) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION PI(49) ,PSI(49) ,NA(2) ,NPSI(2) 

N=NA( 1) 

L=1 


OP 03010 
OP 03020 
OP 03030 
OP 03040 
OP 03050 
OP 03060 
OP 03070 
OP 03080 
OP 03090 
OP 03100 
OP 03110 
OP 03120 
OP 03130 
OP 03140 
OP 03150 
OP 03160 
OP 03170 
OP 03180 
OP 03190 
OP 03200 
OP 03210 
OP 03220 
OP 03230 
OP 03240 
OP 03250 
OP 03260 
OP 03270 
OP 03280 
OP 03290 
OP 03300 
OP 03310 
OP 03320 
OP 03330 
OP 03340 
OP 03350 
OP 03360 
OP 03370 
OP 03380 
OP 03390 
OP 03400 
OP 03410 
OP 03420 
OP 03430 
OP 03440 
OP 03450 
OP 03460 
OP 03470 
OP 03480 
OP 03490 
OP 03500 
OP 03510 
OP 03520 
OP 03530 
OP 03540 
OP 03550 
OP 03560 
OP 03570 
OP 03580 
OP 03590 
OP 03600 
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DO 10 1=1, N 

OP 03610 

DO 20 J=1 ,N 

OP 03620 

K=N*(N-I)+J 

OP 03630 

PSI(L)=PI(K) 

OP 03640 

L=L+1 

OP 03650 

20 CONTINUE 

OP 03660 

10 CONTINUE 

OP 03670 

RETURN 

OP 03680 

END 

OP 03690 

C **** SUBROUTINE NORMAL 

OP 03700 

SUBROUTINE NORMAL(A,NA,B,NB) 

OP 03710 

IMPLICIT REAL*8 (A-H.O-Z) 

OP 03720 

DIMENSION A(49) ,B(49) ,C(7) ,NA(2) ,NB(2) 

OP 03730 

C COMPUTES EUCLIDIAN NORM OF EACH COLUMN 

OP 03740 

N=NA( 1) 

OP 03750 

K=0 

OP 03760 

DO 10 1=1, N 

OP 03770 

SUM=0. 

OP 03780 

DO 20 J=1,N 

OP 03790 

J1=J+K 

OP 03800 

TEMP=A(J1)*A(J1) 

OP 03810 

SUM=SUM+TEMP 

OP 03820 

20 CONTINUE 

OP 03830 

K=K+N 

OP 03840 

C(I)=DSQRT(SUM) 

OP 03850 

10 CONTINUE 

OP 03860 

C NORMALIZE EACH COLUMN 

OP 03870 

K=0 

OP 03880 

DO 30 1=1, N 

OP 03890 

DO 40 J=1,N 

OP 03900 

J1=J+K 

OP 03910 

B(J1)=A(J1)/C(I) 

OP 03920 

40 CONTINUE 

OP 03930 

K=K+N 

OP 03940 

30 CONTINUE 

OP 03950 

RETURN 

OP 03960 

END 

OP 03970 



